Arcsine Distribution#
The arcsine distribution is a continuous distribution on a finite interval whose density spikes near the endpoints.
It shows up in two classic ways:
As a special Beta distribution: (\text{Beta}(\tfrac12,\tfrac12)) on ([0,1]).
As the limit law in the arcsine laws for symmetric random walks / Brownian motion (e.g., the fraction of time spent above 0).
Learning goals#
By the end, you should be able to:
write down and interpret the pdf/cdf of the arcsine distribution
derive its mean and variance using a clean trig substitution
sample it efficiently with NumPy only (inverse CDF)
use
scipy.stats.arcsineforpdf,cdf,rvs, andfitrecognize where it appears in statistics (Jeffreys prior) and stochastic processes (arcsine laws)
Prerequisites#
Calculus: change of variables, basic integrals
Probability: pdf/cdf, expectation and variance
Familiarity with the Beta distribution is helpful but not required
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from scipy import special, stats
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(42)
1) Title & Classification#
Distribution name: arcsine
Type: continuous
Support:
Canonical form: (x \in [0,1])
General (endpoint) form: (x \in [a,b]) with (a<b)
Parameter space:
Endpoint form: ((a,b) \in \mathbb{R}^2) with (a<b)
SciPy form (location/scale):
loc = a,scale = b-awithscale > 0
Key identity:
[ X \sim \text{arcsine on }[0,1] \quad\Longleftrightarrow\quad X \sim \text{Beta}\left(\tfrac12,\tfrac12\right). ]
2) Intuition & Motivation#
What it models#
The arcsine distribution is a natural model for fractions and time proportions that tend to be near 0 or near 1 more often than “in the middle”.
Its pdf is U-shaped on ([0,1]):
very high density near 0 and 1
lowest density near 1/2
This makes it a good model when boundary behavior is common (e.g., “almost always” vs “almost never”).
Where it comes from (a memorable construction)#
Let (\Theta \sim \mathrm{Unif}(0,\pi)). Then
[ Y = \cos\Theta \in [-1,1] ]
has density (f_Y(y)=\tfrac{1}{\pi\sqrt{1-y^2}}), the classic arcsine law on ([-1,1]).
If we map ([-1,1]) to ([0,1]) via (X = \tfrac{Y+1}{2}), we obtain
[ f_X(x) = \frac{1}{\pi\sqrt{x(1-x)}}, \qquad x\in(0,1). ]
This trig construction is also the cleanest route to the mean/variance derivations later.
Typical real-world use cases#
Stochastic processes (arcsine laws): for a symmetric random walk / Brownian motion, the fraction of time the process stays above 0 has an arcsine limit law.
Bayesian statistics: (\text{Beta}(\tfrac12,\tfrac12)) is the Jeffreys prior for a Bernoulli/Binomial probability parameter (p). It is “noninformative” in an information-geometric sense and is heavier near 0 and 1 than the uniform prior.
Generative modeling of probabilities: when you want random probabilities that are often extreme (near 0 or 1), arcsine is a simple choice.
Relations to other distributions#
Special case of the Beta family: (\text{arcsine} = \text{Beta}(\tfrac12,\tfrac12))
Location/scale transform to any ([a,b])
Closely related to the “arcsine on ([-1,1])” density (\tfrac{1}{\pi\sqrt{1-y^2}})
3) Formal Definition#
We’ll use the endpoint parameterization ((a,b)) with (a<b).
PDF#
[ f(x\mid a,b) = \frac{1}{\pi\sqrt{(x-a)(b-x)}}, \qquad x\in(a,b) ]
and (f(x\mid a,b)=0) for (x\notin[a,b]).
Notes:
The pdf diverges as (x\to a) or (x\to b), but the divergence is integrable.
In SciPy, this corresponds to
stats.arcsine(loc=a, scale=b-a).
CDF#
[ F(x\mid a,b) = \begin{cases} 0, & x\le a,\ \frac{2}{\pi}\arcsin!\Big(\sqrt{\frac{x-a}{b-a}}\Big), & a<x<b,\ 1, & x\ge b. \end{cases} ]
In the canonical ([0,1]) case:
[ f(x)=\frac{1}{\pi\sqrt{x(1-x)}},\quad F(x)=\frac{2}{\pi}\arcsin(\sqrt{x}). ]
def _check_ab(a: float, b: float) -> None:
if not (np.isfinite(a) and np.isfinite(b) and a < b):
raise ValueError("Require finite endpoints with a < b.")
def arcsine_pdf(x, a: float = 0.0, b: float = 1.0):
'''Arcsine pdf on [a,b]. Vectorized over x.'''
_check_ab(a, b)
x = np.asarray(x, dtype=float)
out = np.zeros_like(x, dtype=float)
interior = (x > a) & (x < b)
boundary = (x == a) | (x == b)
out[boundary] = np.inf
out[interior] = 1.0 / (np.pi * np.sqrt((x[interior] - a) * (b - x[interior])))
return out
def arcsine_logpdf(x, a: float = 0.0, b: float = 1.0):
'''Numerically friendlier log-pdf (still +inf at boundaries).'''
_check_ab(a, b)
x = np.asarray(x, dtype=float)
out = np.full_like(x, -np.inf, dtype=float)
interior = (x > a) & (x < b)
boundary = (x == a) | (x == b)
out[boundary] = np.inf
out[interior] = (
-np.log(np.pi)
- 0.5 * np.log(x[interior] - a)
- 0.5 * np.log(b - x[interior])
)
return out
def arcsine_cdf(x, a: float = 0.0, b: float = 1.0):
'''Arcsine CDF on [a,b].'''
_check_ab(a, b)
x = np.asarray(x, dtype=float)
z = (x - a) / (b - a)
z = np.clip(z, 0.0, 1.0) # protects sqrt/arcsin from tiny floating-point drift
out = (2.0 / np.pi) * np.arcsin(np.sqrt(z))
out = np.where(x <= a, 0.0, out)
out = np.where(x >= b, 1.0, out)
return out
def arcsine_ppf(u, a: float = 0.0, b: float = 1.0):
'''Inverse CDF (percent point function).'''
_check_ab(a, b)
u = np.asarray(u, dtype=float)
if np.any((u < 0.0) | (u > 1.0)):
raise ValueError("u must be in [0,1].")
return a + (b - a) * np.sin(0.5 * np.pi * u) ** 2
def arcsine_rvs(size=None, a: float = 0.0, b: float = 1.0, rng: np.random.Generator | None = None):
'''NumPy-only sampling via inverse transform.'''
_check_ab(a, b)
if rng is None:
rng = np.random.default_rng()
u = rng.random(size=size)
return arcsine_ppf(u, a=a, b=b)
# Quick sanity check against SciPy on [0,1]
x_grid = np.linspace(1e-6, 1 - 1e-6, 5)
np.c_[x_grid, arcsine_pdf(x_grid), stats.arcsine.pdf(x_grid)]
array([[ 0.000001, 318.310045, 318.310045],
[ 0.25 , 0.735105, 0.735105],
[ 0.5 , 0.63662 , 0.63662 ],
[ 0.75 , 0.735105, 0.735105],
[ 0.999999, 318.310045, 318.310045]])
4) Moments & Properties#
Let (X \sim \mathrm{Arcsine}(a,b)).
Mean and variance#
[ \mathbb{E}[X] = \frac{a+b}{2}, \qquad \mathrm{Var}(X) = \frac{(b-a)^2}{8}. ]
Because the distribution is symmetric around (\tfrac{a+b}{2}), the skewness is 0.
Kurtosis#
For the canonical ([0,1]) case (equivalently (\text{Beta}(\tfrac12,\tfrac12))):
skewness = 0
excess kurtosis = (-\tfrac{3}{2}) (so kurtosis = (3 - \tfrac{3}{2} = \tfrac{3}{2}))
These values stay the same under location/scale transforms.
MGF and characteristic function#
It’s convenient to express these using Bessel functions. Define:
(I_0): modified Bessel function of the first kind (order 0)
(J_0): Bessel function of the first kind (order 0)
Then
[ M_X(t) = \mathbb{E}[e^{tX}] = \exp\Big( t,\tfrac{a+b}{2}\Big), I_0\Big( t,\tfrac{b-a}{2}\Big), \qquad t\in\mathbb{R} ]
[ \varphi_X(t) = \mathbb{E}[e^{itX}] = \exp\Big( it,\tfrac{a+b}{2}\Big), J_0\Big( t,\tfrac{b-a}{2}\Big), \qquad t\in\mathbb{R}. ]
Entropy (differential)#
For the canonical ([0,1]) distribution:
[ H(X) = \log\Big(\frac{\pi}{4}\Big). ]
On ([a,b]), scaling adds (\log(b-a)):
[ H(X) = \log(b-a) + \log\Big(\frac{\pi}{4}\Big) = \log\Big( (b-a),\frac{\pi}{4}\Big). ]
def arcsine_mean(a: float = 0.0, b: float = 1.0) -> float:
_check_ab(a, b)
return 0.5 * (a + b)
def arcsine_var(a: float = 0.0, b: float = 1.0) -> float:
_check_ab(a, b)
return (b - a) ** 2 / 8.0
def arcsine_mgf(t, a: float = 0.0, b: float = 1.0):
'''MGF using modified Bessel I0.'''
_check_ab(a, b)
t = np.asarray(t, dtype=float)
return np.exp(t * 0.5 * (a + b)) * special.i0(t * 0.5 * (b - a))
def arcsine_cf(t, a: float = 0.0, b: float = 1.0):
'''Characteristic function using Bessel J0.'''
_check_ab(a, b)
t = np.asarray(t, dtype=float)
return np.exp(1j * t * 0.5 * (a + b)) * special.j0(t * 0.5 * (b - a))
def arcsine_entropy(a: float = 0.0, b: float = 1.0) -> float:
_check_ab(a, b)
return float(np.log((b - a) * np.pi / 4.0))
# Monte Carlo check on [0,1]
n_mc = 200_000
x_mc = arcsine_rvs(n_mc, rng=rng)
print("Monte Carlo mean:", x_mc.mean(), "| theory:", arcsine_mean())
print("Monte Carlo var :", x_mc.var(), "| theory:", arcsine_var())
print("Entropy theory :", arcsine_entropy())
Monte Carlo mean: 0.4999078439779247 | theory: 0.5
Monte Carlo var : 0.12498062340322307 | theory: 0.125
Entropy theory : -0.2415644752704905
5) Parameter Interpretation#
The parameters (a) and (b) are endpoints of the support.
(a) shifts the distribution left/right
(b-a) scales (stretches) the interval
Crucially, the shape in normalized coordinates does not change.
If (X \sim \mathrm{Arcsine}(a,b)) and
[ Z = \frac{X-a}{b-a}, ]
then (Z \sim \mathrm{Arcsine}(0,1)) (equivalently (\text{Beta}(\tfrac12,\tfrac12))).
So the family is “rigid”: parameters only translate and rescale.
intervals = [(0.0, 1.0), (-1.0, 1.0), (2.0, 5.0)]
fig = go.Figure()
for a, b in intervals:
x = np.linspace(a + 1e-4 * (b - a), b - 1e-4 * (b - a), 800)
fig.add_trace(
go.Scatter(
x=x,
y=arcsine_pdf(x, a=a, b=b),
name=f"pdf on [{a:g}, {b:g}]",
)
)
fig.update_layout(
title="Arcsine pdf for different endpoints (note the endpoint spikes)",
xaxis_title="x",
yaxis_title="density",
)
fig.show()
6) Derivations#
6.1 Expectation and variance (trig substitution)#
A standard trick is to parameterize (x\in[a,b]) by an angle.
Let
[ x(\theta) = \frac{a+b}{2} + \frac{b-a}{2}\cos\theta,\qquad \theta\in(0,\pi). ]
Then
[ (x-a)(b-x) = \Big(\tfrac{b-a}{2}\Big)^2 \sin^2\theta, \qquad \mathrm{d}x = -\tfrac{b-a}{2}\sin\theta,\mathrm{d}\theta. ]
Plugging into (f(x\mid a,b),\mathrm{d}x):
[ \frac{1}{\pi\sqrt{(x-a)(b-x)}},\mathrm{d}x = \frac{1}{\pi,(\tfrac{b-a}{2})\sin\theta},\Big(-\tfrac{b-a}{2}\sin\theta,\mathrm{d}\theta\Big) = -\frac{1}{\pi},\mathrm{d}\theta. ]
Flipping integration limits turns the minus sign into a plus, so effectively
[ f(x),\mathrm{d}x = \frac{1}{\pi},\mathrm{d}\theta, \quad \theta\sim\mathrm{Unif}(0,\pi). ]
Mean.
[ \mathbb{E}[X] = \frac{1}{\pi}\int_0^\pi x(\theta),\mathrm{d}\theta = \frac{1}{\pi}\int_0^\pi \Big(\tfrac{a+b}{2} + \tfrac{b-a}{2}\cos\theta\Big),\mathrm{d}\theta = \frac{a+b}{2} ]
because (\int_0^\pi \cos\theta,\mathrm{d}\theta = 0).
Variance. Write (\mu=\tfrac{a+b}{2}). Then
[ X-\mu = \tfrac{b-a}{2}\cos\theta, \qquad (X-\mu)^2 = \Big(\tfrac{b-a}{2}\Big)^2 \cos^2\theta. ]
Therefore
[ \mathrm{Var}(X) = \mathbb{E}[(X-\mu)^2] = \Big(\tfrac{b-a}{2}\Big)^2 \cdot \frac{1}{\pi}\int_0^\pi \cos^2\theta,\mathrm{d}\theta = \Big(\tfrac{b-a}{2}\Big)^2 \cdot \frac{1}{2} = \frac{(b-a)^2}{8}. ]
6.2 Likelihood#
For i.i.d. data (x_1,\dots,x_n) with all points in ((a,b)), the likelihood is
[ L(a,b) = \prod_{i=1}^n \frac{1}{\pi\sqrt{(x_i-a)(b-x_i)}}. ]
The log-likelihood (up to an additive constant) is
[ \ell(a,b) = -\frac{1}{2}\sum_{i=1}^n \log(x_i-a) - \frac{1}{2}\sum_{i=1}^n \log(b-x_i), \qquad \text{subject to } a<\min_i x_i,; b>\max_i x_i. ]
Important: because (\log(x_i-a)\to -\infty) as (a\uparrow \min x_i), the log-likelihood can grow without bound. So the unconstrained MLE for ((a,b)) does not exist (the likelihood is unbounded). Practical fitting methods therefore use constraints/regularization or alternative estimators.
7) Sampling & Simulation#
Inverse transform sampling (NumPy-only)#
Starting from the CDF for ([a,b]):
[ u = F(x\mid a,b) = \frac{2}{\pi}\arcsin!\Big(\sqrt{\frac{x-a}{b-a}}\Big) ]
Solve for (x):
[ \sqrt{\frac{x-a}{b-a}} = \sin\Big(\frac{\pi u}{2}\Big) \quad\Rightarrow\quad x = a + (b-a),\sin^2\Big(\frac{\pi u}{2}\Big). ]
Algorithm
Sample (u \sim \mathrm{Unif}(0,1))
Return (x = a + (b-a)\sin^2(\pi u/2))
This is exact, fast, and requires only NumPy.
# Sampling demo
a, b = 0.0, 1.0
x = arcsine_rvs(10, a=a, b=b, rng=rng)
x
array([0.35102 , 0.333278, 0.957626, 0.620713, 0.864426, 0.331615,
0.148247, 0.814071, 0.448362, 0.69346 ])
8) Visualization#
We’ll visualize:
the pdf
the cdf
Monte Carlo samples vs the theoretical pdf
# PDF and CDF on [0,1]
eps = 1e-4
x = np.linspace(eps, 1 - eps, 1000)
fig_pdf = go.Figure()
fig_pdf.add_trace(go.Scatter(x=x, y=arcsine_pdf(x), name="arcsine pdf"))
fig_pdf.add_trace(go.Scatter(x=x, y=np.ones_like(x), name="uniform(0,1) pdf", line=dict(dash="dash")))
fig_pdf.update_layout(title="PDF on [0,1] (arcsine vs uniform)", xaxis_title="x", yaxis_title="density")
fig_pdf.show()
fig_cdf = go.Figure()
fig_cdf.add_trace(go.Scatter(x=x, y=arcsine_cdf(x), name="arcsine cdf"))
fig_cdf.add_trace(go.Scatter(x=x, y=x, name="uniform(0,1) cdf", line=dict(dash="dash")))
fig_cdf.update_layout(title="CDF on [0,1] (arcsine vs uniform)", xaxis_title="x", yaxis_title="F(x)")
fig_cdf.show()
# Monte Carlo samples vs pdf
n = 60_000
samples = arcsine_rvs(n, rng=rng)
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=samples,
nbinsx=80,
histnorm="probability density",
name="samples (hist)",
opacity=0.6,
)
)
fig.add_trace(go.Scatter(x=x, y=arcsine_pdf(x), name="theoretical pdf", line=dict(color="black")))
fig.update_layout(
title="Monte Carlo histogram with theoretical pdf overlay",
xaxis_title="x",
yaxis_title="density",
barmode="overlay",
)
fig.show()
print("sample mean/var:", samples.mean(), samples.var())
sample mean/var: 0.5012310506354244 0.12544470668500474
9) SciPy Integration#
SciPy provides the distribution as scipy.stats.arcsine.
stats.arcsine.pdf(x, loc=a, scale=b-a)stats.arcsine.cdf(x, loc=a, scale=b-a)stats.arcsine.rvs(size=..., loc=a, scale=b-a, random_state=...)stats.arcsine.fit(data)estimateslocandscale
We’ll also verify the identity with (\text{Beta}(\tfrac12,\tfrac12)).
# SciPy: pdf/cdf/rvs
x = np.linspace(1e-4, 1 - 1e-4, 1000)
pdf_scipy = stats.arcsine.pdf(x)
cdf_scipy = stats.arcsine.cdf(x)
# Identity: arcsine == Beta(1/2, 1/2)
pdf_beta = stats.beta(a=0.5, b=0.5).pdf(x)
print("max |pdf_scipy - pdf_beta|:", np.max(np.abs(pdf_scipy - pdf_beta)))
# Sampling
s_scipy = stats.arcsine.rvs(size=5, random_state=rng)
s_numpy = arcsine_rvs(5, rng=rng)
print("SciPy rvs:", s_scipy)
print("NumPy rvs:", s_numpy)
# Fitting loc/scale (note: likelihood is tricky near the endpoints)
data = arcsine_rvs(2_000, a=2.0, b=5.0, rng=rng)
loc_hat, scale_hat = stats.arcsine.fit(data)
print("fit loc, scale:", loc_hat, scale_hat)
print("true loc, scale:", 2.0, 3.0)
max |pdf_scipy - pdf_beta|: 3.552713678800501e-15
SciPy rvs: [0.425294 0.712055 0.864634 0.234072 0.38168 ]
NumPy rvs: [0.638791 0.463073 0.730777 0.226557 0.552068]
fit loc, scale: 2.0000041654452665 3.00059414578592
true loc, scale: 2.0 3.0
10) Statistical Use Cases#
10.1 Hypothesis testing (arcsine law for random walks)#
For a symmetric random walk (S_t = \sum_{i=1}^t \epsilon_i) with (\epsilon_i\in{-1,+1}) i.i.d., one version of the arcsine law says:
the fraction of time the walk is positive converges in distribution to an arcsine law.
We’ll simulate many random walks, compute the proportion of steps with (S_t>0), and compare to the arcsine distribution on ([0,1]).
n_paths = 4000
n_steps = 600
steps = rng.choice([-1, 1], size=(n_paths, n_steps))
paths = np.cumsum(steps, axis=1)
frac_positive = (paths > 0).mean(axis=1)
ks = stats.kstest(frac_positive, stats.arcsine.cdf)
print("KS statistic:", ks.statistic)
print("KS p-value :", ks.pvalue)
x = np.linspace(1e-4, 1 - 1e-4, 800)
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=frac_positive,
nbinsx=60,
histnorm="probability density",
name="random-walk fractions",
opacity=0.6,
)
)
fig.add_trace(go.Scatter(x=x, y=stats.arcsine.pdf(x), name="arcsine pdf", line=dict(color="black")))
fig.update_layout(
title="Fraction of time positive in a symmetric random walk (simulation)",
xaxis_title="fraction of steps with S_t > 0",
yaxis_title="density",
barmode="overlay",
)
fig.show()
KS statistic: 0.02575
KS p-value : 0.009764424160149908
10.2 Bayesian modeling (Jeffreys prior for Bernoulli/Binomial)#
For a Bernoulli/Binomial success probability (p), the Jeffreys prior is
[ p \sim \mathrm{Beta}(\tfrac12, \tfrac12), ]
which is exactly the arcsine distribution on ([0,1]).
If we observe (k) successes in (n) trials, the posterior is
[ p \mid k \sim \mathrm{Beta}\Big(k+\tfrac12,; n-k+\tfrac12\Big). ]
We’ll compare Jeffreys’ prior to a uniform prior (\mathrm{Beta}(1,1)).
n, k = 20, 3
prior_jeffreys = stats.beta(0.5, 0.5)
prior_uniform = stats.beta(1.0, 1.0)
post_jeffreys = stats.beta(k + 0.5, n - k + 0.5)
post_uniform = stats.beta(k + 1.0, n - k + 1.0)
x = np.linspace(1e-4, 1 - 1e-4, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=prior_jeffreys.pdf(x), name="Jeffreys prior Beta(1/2,1/2)", line=dict(dash="dash")))
fig.add_trace(go.Scatter(x=x, y=prior_uniform.pdf(x), name="Uniform prior Beta(1,1)", line=dict(dash="dash")))
fig.add_trace(go.Scatter(x=x, y=post_jeffreys.pdf(x), name=f"Posterior (Jeffreys), k={k}, n={n}"))
fig.add_trace(go.Scatter(x=x, y=post_uniform.pdf(x), name=f"Posterior (uniform), k={k}, n={n}"))
ci_low, ci_high = post_jeffreys.ppf([0.025, 0.975])
fig.add_vline(x=ci_low, line=dict(color="black", dash="dot"))
fig.add_vline(x=ci_high, line=dict(color="black", dash="dot"))
fig.update_layout(
title="Jeffreys (arcsine) prior and resulting posterior for a Binomial proportion",
xaxis_title="p",
yaxis_title="density",
)
fig.show()
print("Jeffreys posterior mean:", post_jeffreys.mean())
print("Jeffreys 95% credible interval:", (ci_low, ci_high))
Jeffreys posterior mean: 0.16666666666666666
Jeffreys 95% credible interval: (0.04413134197515587, 0.348577710858454)
10.3 Generative modeling (sampling extreme probabilities)#
If you sample (p) from an arcsine distribution and then sample data conditional on (p), you generate datasets where the latent probability is often near 0 or 1.
One simple example:
sample (p \sim \mathrm{Beta}(\tfrac12,\tfrac12))
then sample a count (K \mid p \sim \mathrm{Binomial}(n,p))
Compared to a uniform prior over (p), this produces more extreme counts (very small or very large (K)).
n = 50
m = 30_000
p_arcsine = stats.arcsine.rvs(size=m, random_state=rng)
p_uniform = rng.random(size=m)
k_arcsine = rng.binomial(n, p_arcsine)
k_uniform = rng.binomial(n, p_uniform)
fig = go.Figure()
fig.add_trace(go.Histogram(x=k_arcsine, histnorm="probability", name="p ~ arcsine", opacity=0.6, nbinsx=n + 1))
fig.add_trace(go.Histogram(x=k_uniform, histnorm="probability", name="p ~ uniform", opacity=0.6, nbinsx=n + 1))
fig.update_layout(
title=f"Counts from Binomial(n={n}, p) with different priors over p",
xaxis_title="k (number of successes)",
yaxis_title="probability",
barmode="overlay",
)
fig.show()
print("P(k=0) arcsine vs uniform:", np.mean(k_arcsine == 0), np.mean(k_uniform == 0))
print("P(k=n) arcsine vs uniform:", np.mean(k_arcsine == n), np.mean(k_uniform == n))
P(k=0) arcsine vs uniform: 0.08146666666666667 0.0201
P(k=n) arcsine vs uniform: 0.07906666666666666 0.020866666666666665
11) Pitfalls#
Invalid parameters: always ensure (a<b) (or
scale>0).Endpoint singularities: the pdf diverges at (a) and (b). This is mathematically fine, but numerically:
avoid evaluating the pdf exactly at the endpoints when plotting
prefer
logpdffor likelihood computations
Floating-point drift: for the CDF, expressions like ((x-a)/(b-a)) can be slightly below 0 or above 1 due to rounding; clipping prevents
sqrt/arcsinfrom producing NaNs.Fitting endpoints: the likelihood can become unbounded as endpoints approach the sample min/max, so naive MLE fitting is ill-posed without additional constraints.
# Numerical illustration: pdf spikes and logpdf is usually safer
x_test = np.array([0.0, 1e-12, 0.5, 1 - 1e-12, 1.0])
print("x:", x_test)
print("pdf:", stats.arcsine.pdf(x_test))
print("logpdf:", stats.arcsine.logpdf(x_test))
x: [0. 0. 0.5 1. 1. ]
pdf: [ inf 318309.886184 0.63662 318313.407023 inf]
logpdf: [ inf 12.670781 -0.451583 12.670792 inf]
12) Summary#
The arcsine distribution on ([0,1]) is (\mathrm{Beta}(\tfrac12,\tfrac12)) with pdf (f(x)=\tfrac{1}{\pi\sqrt{x(1-x)}}).
It has a U-shaped density: lots of mass near the boundaries.
Endpoint form on ([a,b]): (f(x\mid a,b)=\tfrac{1}{\pi\sqrt{(x-a)(b-x)}}).
Mean (=\tfrac{a+b}{2}), variance (=\tfrac{(b-a)^2}{8}); skewness 0, excess kurtosis (-\tfrac{3}{2}).
Fast exact sampling: (x = a + (b-a)\sin^2(\pi u/2)) with (u\sim\mathrm{Unif}(0,1)).
In practice,
scipy.stats.arcsinegivespdf,cdf,rvs, andfit(with caution near endpoints).